#ifndef BTL_MATHS_LAPACKSVD
#define BTL_MATHS_LAPACKSVD

#include <Eigen/Core>

namespace btl
{
namespace maths
{

class LapackSVD {
   public:
      typedef double Scalar;

      LapackSVD();

      template <typename MT>
      LapackSVD(const Eigen::MatrixBase<MT>& m, int opts = 0);

      template <typename MT>
      void compute(const Eigen::MatrixBase<MT>& m, int opts = 0);

      const Eigen::VectorXd& singularValues() const;
      const Eigen::MatrixXd& matrixU() const;
      Eigen::MatrixXd::ConstTransposeReturnType matrixV() const;

   private:
      void doCompute(Eigen::MatrixXd& mat, int opts);

      Eigen::VectorXd _sv;
      Eigen::MatrixXd _VT;
      Eigen::MatrixXd _U;
      int _opts;
};

} // namespace maths
} // namespace btl

namespace btl
{

using maths::LapackSVD;

} // namespace btl

// ====================================================================
// === Implementation

namespace btl
{
namespace maths
{

inline LapackSVD::LapackSVD()
{
}

template <typename MT>
inline LapackSVD::LapackSVD(const Eigen::MatrixBase<MT>& m, int opts)
{
   compute(m, opts);
}

template <typename MT>
inline void LapackSVD::compute(const Eigen::MatrixBase<MT>& mat, int opts)
{
   Eigen::MatrixXd mfull = mat;
   doCompute(mfull, opts);
}

inline const Eigen::VectorXd& LapackSVD::singularValues() const
{
   return _sv;
}

inline const Eigen::MatrixXd& LapackSVD::matrixU() const
{
   assert(_opts & (Eigen::ComputeThinU | Eigen::ComputeFullU));
   return _U;
}

inline Eigen::MatrixXd::ConstTransposeReturnType LapackSVD::matrixV() const
{
   assert(_opts & (Eigen::ComputeThinV | Eigen::ComputeFullV));
   return _VT.transpose();
}

} // namespace maths
} // namespace btl

#endif // BTL_MATHS_LAPACKSVD
